Real-time wavelet-transform spectrum analyzer for the investigation of 1//" noise 
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A wavelet transform spectrum analyzer operating in real time within the frequency range 3 x 
10~^ — 1.3 X lO^Hz has been implemented on a low-cost Digital Signal Processing board operating 
at 150MHz. The wavelet decomposition of the signal allows to efficiently process non-stationary 
signals dominated by large amplitude events fairly well localized in time, thus providing the natural 
tool to analyze processes characterized by 1//" power spectrum. The parallel architecture of the 
DSP allows the real-time processing of the wavelet transform of the signal sampled at 0.3MHz. The 
bandwidth is about 220dB, almost ten decades. The power spectrum of the scattered intensity is 
processed in real time from the mean square value of the wavelet coefficients within each frequency 
band. The performances of the spectrum analyzer have been investigated by performing Dynamic 
Light Scattering experiments on colloidal suspensions and by comparing the measured spectra with 
the correlation functions data obtained with a traditional multi tau correlator. In order to asses the 
potentialities of the spectrum analyzer in the investigation of processes involving a wide range of 
timescales, we have performed measurements on a model system where fluctuations in the scattered 
intensities are generated by the number fluctuations in a dilute colloidal suspension illuminated by 
a wide beam. This system is characterized by a power-law spectrum with exponent —3/2 in the 
scattered intensity fluctuations. The spectrum analyzer allows to recover the power spectrum with 
a dynamic range spanning about 8 decades. The advantages of wavelet analysis versus correlation 
analysis in the investigation of processes characterized by a wide distribution of time scales and 
non-stationary processes are briefly discussed. 

PACS numbers: 07.50.Qx, 07.60. Rd, 02.70.Rr, 05.40.Ca, 05.45.Df 



I. INTRODUCTION 

The investigation of random processes characterized 
by a wide range of timescales is becoming increasingly 
important in many experimental fields ranging from the 
investigation of earthquakes in earth-science, to turbu- 
lence in fluids, or the price fluctuations in flnancial mar- 
kets H, ^, ^. In particular, many random processes 
exhibit a scale invariant structure. The self-similarity of 
the signal is reflected in the power law behavior of the 
power spectrum, which lacks characteristic time scales. 

Traditional spectral decomposition techniques often 
fail when applied to a self similar signal. This is due 
to the fact that such a signal is typically the superpo- 
sition of bursts occurring at many different timescales. 
Although rare, long bursts having a very large ampli- 
tude provide much of the energy content of the signal, 
while short bursts with small amplitude, although very 
frequent, give small contribution to the energy of the sig- 
nal. The energy is therefore fairly well localized in time 
and the signal is non-stationary. This localization pre- 
vents Fourier Analisys to work as an effective tool for 
this kind of signal. 

However, a scale invariant signal can be efficiently an- 
alyzed by using the renormalization group approach de- 
vised in statistical physics to describe systems at a critical 
point . The self similar system goes through a consec- 
utive sequence of coarse grainings which generate a new 
signal with statistical properties similar to those of the 
original one. 

In practice the renormalization group analisys of a sig- 



nal is achieved by using the wavelet transform technique. 
Basically the signal is convolved with a chosen mother 
wavelet function. The mother wavelet consists in a trun- 
cated wave packet with central frequency w and band- 
width r, localized at time t, with length / cx l/F. The 
mean square amplitude of this component provides the 
power spectral amplitude of the signal in the frequency 
range (w — F, a; -I- F). An in-quadrature component of 
the signal is also simultaneously processed to recover 
the discarded lower frequency component, representing 
the coarse grained signal. The mother wavelet is then 
rescaled on a larger scale and the process iterated. 

In this paper we describe a real-time wavelet transform 
spectrum analyzer (WTSA) working in the frequency 
range 3 x 10~^ — 1.3 x lO^Hz. The spectrum analyzer 
is built around a Texas Instruments C6711 DSP Starter 
Kit (DSK). This board incorporates a C6711 DSP run- 
ning at ISOMHz together with a parallel port to commu- 
nicate with a Personal computer and various connectors 
to interface the board to the real world. Solutions based 
on wavelet analysis by means of the post processing on a 
dedicated personal computer of the digitized signal have 
recently been proposed |^. The highly parallel pro- 
cessing architecture of the DSP allows the real-time pro- 
cessing of the wavelet transform of the signal sampled at 
O.SMHz. The bandwidth spans almost ten decades. 

We will present applications to Dynamic Light Scatter- 
ing (DLS) from colloidal suspensions and compare the re- 
sults obtained with the WTSA with the correlation func- 
tion obtained with a traditional multi tau correlator. Dy- 
namic Light Scattering has been used very extensively to 
investigate stochastic fluctuations in the intensity of the 
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light scattered by suspended particles undergoing brow- 
nian motion, aggregation processes, gelation processes, 
critical dynamics ^ |l^. In recent times, processes 
characterized by a wide distribution of relaxation times, 
often leading to stretched-exponential decay of the au- 
tocorrelation function, have attracted much interest (see 
Refs. ^ and references therein). The development 
of log-scale correlators has allowed the investigation of 
processes extending from about 100ns up to hours. How- 
ever, large time lags require very long acquisitions times 
to get good statistical accuracy in the correlation tails. 

The wavelet analisys of the fluctuations in the intensity 
of the scattered light provides higher accuracy in shorter 
times when compared to DLS when non-stationary self 
similar signals are processed. This is due to the fact that 
the localized wavelets used to process the signal are very 
efficient in isolating the well localized features of a self- 
similar signal, and this provides a shorter convergence 
time of the processing to get results with a chosen statis- 
tical accuracy. 

To asses the potentialities of the WTSA in investigat- 
ing time-scale invariant systems we have devised a simple 
model system giving rise to a power-law spectrum with 
exponent —3/2. This system is made up by strongly 
diluted brownian particles diffusing within a large sam- 
pling volume. We show that the WTSA is able to mea- 
sure the power spectrum within a dynamic range as large 
as 8 decades. We also show that the WTSA requires a 
measurement time roughly 100 times smaller than that 
needed by a correlator to get the same statistical accu- 
racy at large lag times. 



II. WAVELET ANALYSIS 

The wavelet analisys of signals is becoming increas- 
ingly popular within the scientific community. With re- 
spect to traditional Fourier decomposition techniques, 
wavelet analisys allows to characterize signals having a 
large bandwidth by using a small number of coefficients. 

Basically, instead of decomposing the signals into 
armonic oscillations as in Fourier analisys, the signal 
gets decomposed into the superposition of wave packets 
(WP). The WP are localized both in time and in fre- 
quency. These wave packets are obtained by rescaling a 
chosen Mother Wavelet (MW) ^' (t) with a proper scale 
s so as to obtain a wavelet function ipto,s {t) — ^' {^~r^)- 
The wavelets exhibit the same functional form of the 
MW, but are centered around time to and their band- 
width is a factor s larger than that of the MW. Therefore, 
if the signal has a bump at time t with a certain dura- 
tion At the wavelet decomposition will likely give rise to 
a wave packet localized at t with a scale s proportional to 
the duration of the bump. This allows a very efficient de- 
composition of spiky signals and non-stationary signals 
in general. In practice the wavelet analysis is usually 
performed by using the popular Discrete Wavelet Trans- 
form (DWT) algorithm |l^. To apply this algorithm. 



one starts with a sampled discrete signal made up by 
2^ samplings, where N is an integer. The signal under- 
goes successive coarse grainings and band-pass filterings 
by convolution with a chosen wavelet set. The scale s of 
the wavelets used during the coarse grainings changes by 
a factor 2 at each step, so that the scaling of the signal 
at the k-th step is given by Sk = sq2''. On the contrary, 
the number of samplings in the coarse grained signal de- 
creases as 2^/2*-'. Therefore, as the scale increases due to 
the successive scale doublings, the resolution in the tem- 
poral localization of the wavelet decreases accordingly. 
In the following we will give a brief and more rigorous 
overview on how the power spectrum is processed from 
the sampled signal. 

Consider a discrete signal sampled at a given sample 
frequency f^. According to Nyquist theorem the higher 
angular frequency it contains is luq — '^fs, that corre- 
sponds to the function = (—1)". The original signal 

is passed through a high-pass filter, giving the wavelet 
coefficients c° , representing the high frequency details of 
the signal in the range (0^0/2, ujq), and through a low-pass 
filter, giving a coarse grained signal y^. Both the filters 
are implemented by convolving the signal with suitable 
filter coefficients gn (the MW) and /i„; they are localized 
in n so that the only non vanishing elements are the ones 
with < n < Af : 

f c" E„ yn9n-2j ^^ 
I y} = E„ Vnhn~2] 

Both the filtered signals c° and are sampled at half the 
original frequency fg. The two filter coefficients hn and 
gn are such that the coarse grained signal y\ , along with 
the details c°, allow to recover the whole original sig- 
nal y^: they are referred to as quadrature mirror filters. 
In the measurements we performed we used Daubechies 
wavelets p^ . 

The procedure is then iterated, using as the new sig- 
nal; a high frequency component is obtained together 
with a coarse grained signal y"^, and so on. In general, 
the following recursion formulas apply: 

At step k, the signal is coarse grained to obtain 
y^i^^; the high frequency details, spanning the angular 
frequency range (u}q/2''~^^ ,ujq/2'') are the wavelet coeffi- 
cients c^j. 

The power spectrum of the signal can then be derived 
from the mean square average of the details (^|c^|^^ ■ 

With the traditional Fourier transform, the analyzed 
wavelengths are discretized, due to the finiteness of the 
sample; the spacing between two consecutive wavenum- 
bers is constant. Unlikely, with wavelet analysis, the 
spacing is exponential: the mean square amplitude of 
the wavelet coefficients represents the power contained 
in an octave. 
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The exponential spacing is particularly suited to mea- 
sure phenomena spanning many decades in wavenumber. 
This feature is analogous to the multi tau analysis done 
by the best commercial correlators. Though the power 
spectrum and the correlation function are connected by 
the well known Wicncr-Kintchinc theorem, in many ac- 
tual cases it is not possible to recover the power spectrum 
by the measurement of the correlation function. In gen- 
eral, the exponential spaced correlation function given 
by a multi-tau correlator cannot be safely Fourier trans- 
formed into the power spectrum, over many decades. 
Moreover, in some cases Wiener-Kintchine theorem does 
not apply. This is the case of some scale invariant and 
non-stationary signals; in Sect. VB we present an exam- 
ple in which the power spectrum obtained by the wavelet 
analysis correctly shows a power law behaviour, while 
a multi tau correlator seemingly says that the signal is 
delta correlated. 

The analogy between the wavelet analysis and the 
multi tau correlation also includes the possibility of per- 
forming the required processing in real time, as the sam- 
ple is acquired. With Fourier analysis, as the sample 
becomes larger, we must analyze more and more wave- 
lengths, and this requires that all the sample is recorded 
and re-analyzed. On the contrary, as time goes on, the 
multi tau correlator averages the older samples; as a new, 
longer delay time is available, the correlation function is 
evaluated from only a few registers. This allows to ex- 
tend the dynamic range in time delay over many decades, 
without the need to record the samples with the same, 
huge dynamic range. The same thing happens with the 
WTSA. The iterative procedure is implemented by using 
TV, almost identical levels k. The memory needed scales 
as while the time needed by the processing does not 
change at all, as we will show in Sect. IV; on the other 
hand, the frequency range scales as 2^. 



III. THE WAVELET TRANSFORM SPECTRUM 
ANALYZER 

The spectrum analyzer is built around a Texas Instru- 
ments C6711 DSP Starter Kit (DSK), graciously pro- 
vided by Texas Instruments DSP University Program. 

This board incorporates a C6711 DSP running at 
150MHz, together with a parallel port to communicate 
with a Personal computer and various connectors to in- 
terface the board to the real world. The software can 
be developed on a Personal Computer by using TI Code 
Composer Studio and then downloaded on the board. 
The board and the bundled software are very cheap com- 
pared to a real time correlator board, the purchase price 
being of the order of 300$ for educational institutions. 

The core of the board is made up by a Digital Sig- 
nal Processor, a single-chip programmable device includ- 
ing a CPU specifically suited for digital signal processing 
0. The DSP is a TMS320C6711, built by Texas Instru- 
ments, one of the fastest processors currently available. 



It is based on VelociTI, a high-performance, advanced 
very-long-instruction-word architecture, and its CPU can 
perform operations on double precision (64 bit) floating 
point numbers. 

The CPU is divided into two almost identical sides; 
each of them includes 16 registers, one word (32 bit) 
long. Couples of registers can be used to represent one 
floating point, double precision register. Each section in- 
cludes four ALUs, operating simultaneously; each ALU 
can perform only a given set of operations. Data paths 
from registers to memory and from a register in one side 
to an ALU in the other are also provided. 

The CPU works at 150MHz. At each clock cycle, up to 
eight instructions are fetched, dispatched and decoded, 
that is, up to one instruction for each ALU. All the 
fixed point instructions are executed in one clock cycle; 
this means that the CPU can perform 1.2 x 10^ instruc- 
tions per second, a figure comparable to the number of 
MFLOPS available on a 2GHz Pentium IV processor. 
Some floating point instructions require more than one 
clock cycle to return the results (execution latency) ; gen- 
erally, the same unit can start a new instruction after 
a few clock cycles (functional unit latency), before the 
results of the previous operations are availabe. For ex- 
ample, the functional unit latency of a double precision 
floating point multiplication is 4 clock cycles, and the 
execution latency is 9 clock cycles. Consider, for ex- 
ample, an ALU starting a multiplication at t = 0: at 
t — A clock cycles it can start a new multiplication and 
at i = 9 clock cycles the results of the first multiplication 
are available on the registers. Since there are two ALUs 
able to perform floating point multiplications, the CPU 
can do 60 x 10® multiplications per second, and, simulta- 
neously, much more floating point additions, using other 
two ALUs. This structure requires a sophisticated, par- 
allel programming technique, resulting in extremely fast, 
real time programs. 

The DSP includes a two level memory |]l^. The flrst 
level is composed by two 4Kbytes cache memories, LIP 
for the program instructions and LID for the data. The 
second level, L2, is a 64Kbytes RAM, that can be config- 
ured partially as a SRAM or as a chache for an external 
memory. Since our program and data fit in the SRAM, 
we didn't use the slower external RAM and we config- 
ured all the L2 as SRAM. Reads and writes in the LIP 
and LID memories are performed in four clock cycles; 
obviously, L2 is slower and a miss in the LI cache results 
in a CPU staU. 

The DSP includes many other peripherals; among 
them, two counter-timers, TIMERl and TIMER2 
They are broadly configurable through a configuration 
register. The input can be connected to the system clock 
or to an external pin, driven by a CMOS signal, thus 
selecting if the device is acting as a counter of external 
events or as a timer. The count register can be read and 
used by the program. Based on the counts and the value 
of a period register, an output is driven, in order to gen- 
erate square waves of different duty cycle. Moreover, the 
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output can drive an interrupt. We used TIMERl as a 
counter of external events. In this way, the input of the 
spectrum analyzer can be driven by a discrete pulse train, 
as in the application to DLS discussed below, where the 
timer is connected to the output of a photon counting 
photomultiplier. Alternatively, by using a voltage to fre- 
quency converter, the input of the tinier could be driven 
by an analog signal, such as the one from a photodiode. 
The other timer, TIMER2, generates a periodic inter- 
rupt, in order to define the integration period: the in- 
terrupt service routine reads the counts of TIMERl and 
resets it, passing the number of counts to the processing 
program. 

Texas Instruments provides a DSK (DSP Starter Kit), 
including a board with the DSP and some peripherals; 
among them, a device that allows to connect the board 
to a PC through a parallel port. The PC can read and 
write the memory of the DSP, without halting the CPU, 
and can run a program. The PC uploads the program 
that evaluates the power spectrum of the signal in the 
SRAM, then runs it, and periodically reads the results. 

The software consists of two parts, one for the DSP 
and one for the host computer. Texas Instruments pro- 
vides Code Composer Studio (CCS), that includes a C 
compiler and an assembler for the DSP. The program for 
the DSP has been written in assembler code, since in this 
way a high processing speed is obtained in conjunction 
with an accurate timing. Moreover, the provided soft- 
ware includes some libraries for building host programs 
that interact with the DSP, that is, to upload a program, 
run it, and read some memory locations in the DSP 16| . 
All the software we developed is freely available. jl7 



IV. DISCRETE WAVELET TRANSFORM 
PROGRAM 

With respect to the hardware correlators traditionally 
used within the Dynamic Light Scattering community, 
the use of a DSP allows the flexible implementation of 
many processing algorithms simply by dowloading a dif- 
ferent processing code onto the board. These algorithms 
include FIR and IIR filters, cross correlation, convolution 
filters etc. 

In this section, we discuss the implementation of the 
discrete, real time wavelet algorithm on the DSK board. 

The overall data structure is shown in Fig. 1, along 
with a simplified description of the data flow. The higest 
sampling frequency is about 0.3MHz and the lowest fre- 
quency corresponds to a period of some hours. Every fre- 
quency octave represents a level in the wavelet processing 
algorithm. Figure 1 shows only levels and 1 and part 
of level 2. The entire structure is made up by iV = 32 
levels, each corresponding to a measured frequency. Ev- 
ery level contains a circular queue of double precision 
floating points to store the last M values of the coarse 
grained signal, where M is the length of the MW. The 
N queues are managed by using the circular addressing 



mode of the DSP; this allows only values of M that are 
powers of two. Therefore, the Daubechies wavelets our 
program can use are those with length M — 4, 8, 16, 32 
and 64. We will discuss the implementation of the Daub4 
Af = 4 transform. 

Level is driven by TIMERl, configured as a counter, 
which counts the digitized pulses at the input of the spec- 
trum analyzer within a selected bin. The integration time 
is defined by TIMER2, configured as a timer, so that its 
output signal generates an interrupt at /<, = 375KIIz, 
corresponding to the sampling frequency. At each sam- 
pling time the number of pulses counted by TIMERl is 
fed to the first element of the shift register y° in Level 
and the other cells are shifted down. Every two sample 
times the content of the shift register is convolved with 
the two wavelet coefficient banks gn and /i„, the MW 
and the smoothing filter coefficients respectively. The 
result of the convolution with the bank hn feeds the shift 
register in Level 1 and gives rise to a coarse grained 
replica of the signal at half the sampling frequency. The 
in-quadrature component obtained from the convolution 
with bank gn contains the wavelet coefficients which 
specify the local amplitude of the band-pass filtered sig- 
nal in the frequency range {ujo/2,ujq), where ujq = irfs. 
The process is repeated in cascade for all the levels. 

The processing of N levels requires to store additional 
informations: the sum and the counter needed to eval- 
uate the mean of |cjj| , the current address inside the 
queue and the status of the queue. 

The coefficients feed the register X] l^nl , which 
stores the mean square amplitude of the high-pass fil- 
tered signal representing the power spectrum S (uj) in 
the angular frequency range (wo/2'^"'"^, wq/Z*^). 

The status of the queue can be (the queue has just 
been processed), 1 (the queue had one input) or 2 (the 
queue is full and requires processing); it is incremented 
by one each time an element is put in the queue and is 
reset to zero when it is processed. At the initialization, 
the status is set to —M + 2, so that the queue is filled 
at least with M elements before it is processed the first 
time. 

As outlined above. Level 1 is updated only when two 
new samplings are available in Level 0. Therefore the 
update of Level 1 occurs with a frequency /s/2, and its 
processing with a frequency /s/4, which corresponds to 
the update frequency of Level 2, and in general Level k is 
updated at a frequency /s/2'' and processed at fs/'i,^'^^- 
The sum of all the processing frequencies for every queue, 
X„/s/2''+^, converges to /s as TV — > -l-oo. This means 
that the frequency at which all the queues are processed 
is smaller than the sampling frequency fs- At each sam- 
pling time, our algorithm processes the lowest-level full 
queue, if any. This ensures that, by processing exactly 
one queue at each sampling time, the output values yfj 
never go to an already full queue. For the Daubechies 
wavelets with M = 4 all the operations can be performed 
in about 2/is. 

Beyond the processing sofware on the DSP board two 
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FIG. 1: Data structure and data flow of the wavelet transform spectrum analyzer. 



host programs run on the host computer: the first up- 
loads the DSP programs and runs it; the second reads 
the memory of the DSP, without halting the CPU. This 
program is used to get a real time display of the processed 
spectrum. The read program accesses the DSP memory 
through the parallel port. Inside the DSK board, data are 
read from memory by using a low priority DMA channel, 
and thus without interfering with the CPU operations. 

Another processing method, not yet implemented, in- 
volves the processing of the cross spectrum of two signals. 
This is very useful to get rid of noise sources such as af- 
terpulsing from a photomultiplier detector. This method 
involves the parallel use of two processing structures like 
the one shown in Fig. 1. Each detector signal feeds one of 
the TIMER input. The output of each structure is a set 
of wavelet coefficients cj„ and c^^^ . The averaged product 
of the coefficients J2 '^an'^tn gives the desired cross spec- 
trum. Due to the simultaneous operation of two pro- 
cessing structures, the minimum sampling frequency is 
reduced by a factor of two. 




FIG. 2: Diagram of the experimental setup. 



V. EXPERIMENTAL RESULTS 

To assess the performances of the Spectrum Analyzer 
we have performed test measurements on a colloidal sus- 
pension and on a model system giving rise to a power-law 
power spectrum. 

The experimental setup is based on the traditional 
DLS one. A diagram of the setup is shown in Fig. 2. 
The beam coming from a JDS Uniphase 35mW CW HeNc 
LASER model 1145P is spatially filtered and focused by 
the lens F at the center of the sample contained in the 
cuvette C. Light scattered from the sample is collected 
by the lens L and by the diaphragms DI and D2. The 
scattered radiation eventually impinges onto the photo- 
cathode of a EMI 9863 B04/350 photon counting Pho- 



tomultiplier Tube (PMT). The scattering angle can be 
changed by rotating the collecting optics and the PMT 
around the cuvette by means of a goniometer. The cur- 
rent pulses from the PMT are amplified and filtered by an 
ALV PM-PD preamplifier-discriminator, which also dig- 
itizes the pulse train at TTL levels. In the traditional 
DLS apparatus the signal is then fed to an ALV5000 
multi tau hardware correlator hosted in a personal com- 
puter. Alternatively, the digitized signal at the output 
of the amplifier-discriminator can be routed to a timer 
input of the DSP board, which is connected to a per- 
sonal computer for the configuration and visualization 
of measurements. The DSP board is also connected to 
a HAMEG 1.6 GHz frequency counter model HM8021-3 
and a HAMEG 20 MHz pulse generator model HM8035 
for testing purposes. 
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FIG. 3: Correlation function of the intensity fluctuations 
of light scattered at 30° by a colloidal suspension, in 
the Dynamic Light Scattering regime, measured with ALV 
equipement. Circles represent the experimental data. The 
continuous line represents a fit of the experimental data up 
to the second cumulant. 



A. Colloidal suspensions 

The test on the colloidal suspension were performed 
on calibrated polystyrene latex of 0.115/zm diameter, di- 
luted at a volume fraction of about = 1.5 x 10^^. For a 
dilute suspension of monodisperse brownian particles the 
Homodyne time autocorrelation function of the intensity 
fluctuations in the scattered light has the usual exponen- 
tial form C (t) = {if [1 + Aexp {-Tt)], where T = 2Dq^ 
and D = ksT/ {6Trria) is the Stokes-Einstein Diffusion co- 
efficient, J] is the shear viscosity for the solvent and a the 
particle radius. The correlation function measured at a 
scattering angle of 30° is plotted in Fig. 3 as a function 
of the lag time. The experimental data show the decay of 
the correlation function within a dynamic range of about 
2.5 decades. The run duration was 2700s. 

At the longer delay times the correlation function is 
dominated by noise due to the poor statistical sample 



accumulated. The additional contribution at small lag 
times is due to the afterpulsing of the PMT. The after- 
pulses of the PMT used to performed the measurement 
have a characteristic range of correlation times between 
O.l^s and 3.2/xs, which give rise to a stretched exponen- 
tial decay at small lags. The decay at lags larger than 
about 3.2/is is determined by the intensity fluctuations 
due to the brownian motion of the colloid. Data can 
be fit with an exponential function exp(— Ft). A more 
refined analisys of the experimental data involves the in- 
clusion of polidispersity effects in the fitting procedure. 
This is traditionally performed by a cumulant analisys 
of the correlation function In the case of a small 
polidispersity the correlation function up to the second 
cumulant is given by 
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exp -Fo |<| + -crH' 



exp (— F |i|) exp 



(r-Fo)^ 



20-2 



dF, 



where the second equality shows that polidispersity ef- 
fects give rise to the superposition of exponential decays 
centered around a central linewidth Fq, these decays be- 
ing weighted by a gaussian function with variance a rep- 
resenting the polidispersity. The best fit of the experi- 
mental data with Eq. ^) is indicated by the solid line in 
Fig. 3 and yields Fq 140Hz and a = 42Hz. 

We performed the same experiment, by using the DSP 
WTSA and the same homodyne detection scheme out- 
lined above. According to the Wiener-Kintchine theo- 
rem, the power spectrum of the fluctuations in the inten- 
sity of light scattered by the colloidal suspension is the 
Fourier transform of the correlation function C (t) . For 
diluted monodisperse brownian particle the power spec- 
trum has a Lorenzian shape: S (uj) oc F/ (F^ -|- w^). 

Figure 4 shows data for the measured power spectrum 
as a function of angular frequency uj. Raw data are rep- 
resented with crosses. At the smaller frequency the effect 
of noise due to the reduced statistical sample is evident. 
This occurs at angular frequencies roughly smaller than 
Irad/s, corresponding to a time delay of the order of 
6s. However, this same noise affects the correlation func- 
tion shown in Fig 3 starting from about t — 10~^s. As 
the duration of the experiment is in both cases 2700s, 
the wavelet processing allows a better convergence at 
larger delay times: the timescale range across which the 
power spectrum can be reliably determined is roughly 
two decades more extended than that of the correlator. 
At larger frequency the power spectrum in Fig. 4 falls 
down up to a constant base line, represented by a hor- 
izontal dashed line. This baseline is determined by the 
shot noise. In fact, it is well known that the detection 
of photons is a random process characterized by a Pois- 
son distribution and gives rise to a delta-correlated white 
noise, whose power spectrum corresponds to the average 
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By taking the Fourier transform of Eq. (0) we get: 
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FIG. 4: Power spectrum of the intensity fluctuations of light 
scattered at 30° by a colloidal suspension, in the Dynamic 
Light Scattering regime, measured with WTSA. Crosses rep- 
resent the raw data. The dashed line is the baseline due to the 
shot noise. Circles represent the power spectrum obtained by 
subtracting the baseline from the raw data. The continuous 
line is the best fit of the experimental data up to the second 
cumulant. 



number of photons detected within the sampling time 
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In the frequency space this gives rise to a constant term 
(n), represented by the dashed line in Fig. 4, which adds 
up to the power spectrum. The average number of pho- 
tons (n) can be easily processed in real time and sub- 
tracted. 

The circles in Fig. 4 show experimental data for the 
power spectrum after the subtraction of the shot noise 
contribution. The data span a dynamic range of more 
than 4 decades. Data can be represented by a Lorenzian 
curve r/ (w^-i-F^). In this case also, a more refined 
analisys involves the inclusion of polidispersity effects. 



Therefore the power spectrum is the superposition of 
lorenzians weighted by the same gaussian function in- 
troduced to describe the polidispersity in the correlation 
analysis. We evaluated numerically the integral of Eq. 
(Q), in order to fit the results in Fig. 4. However, each cir- 
cle in Fig. 4 represents the integral of the power spectrum 
across an octave; this smooths even more the bell-shaped 
measured curve; this effect has been considered in the fit 
procedure. We obtained the best fit with Fq = ISSHz 
and a — 37IIz, which favourably compare with the re- 
sults obtained with the correlator. The best fitting curve 
is represented by the solid line. 

At the higher frequencies the spectrum is dominated 
by the contribution of afterpulses to the intensity fiuctu- 
ations which give rise to the tail in the range 10'' — lO^Hz 
in Fig. 4. 

The wavelet analysis of the data presents some ad- 
vantages when compared to correlation techniques. A 
notable advantage is represented by the way background 
subtraction is achieved. In correlation techniques a con- 
stant baseline has to be subtracted from the measured 
correlation function to recover data like those shown in 
Fig. 3. The baseline represents the correlation function 
at large la times, where correlations disappear. The in- 
trinsic contrast of the correlation function, defined by 
the ratio between the amplitude of the completely cor- 
related part at i = and the uncorrelated part at large 
lag times, [C (0) — C (-f-oo)] /C (oo), cannot exceed the 
value of two. Therefore, the signal and the background 
are of the same order of magnitude and the dynamic 
range is strongly limited by the accuracy in the base- 
line. The accuracy in the baseline can be determined 
from the arguments presented in Ref. and summa- 
rized below. Let us call tl the largest decay time in the 
correlation function and T the duration of the run. Sam- 
ples acquired at time lags smaller than are correlated, 
therefore the number of independent samples acquired 
during T is = T/tl- This independent samples av- 
erage out according to a normal distribution to give rise 
to the baseline. Therefore, the accuracy in the base- 

1 /2 

line is of the order of 1/iVc . In the case of the results 
presented in Fig. 3, tl = 40ms, T — 2700s, thus giv- 
ing X/nI/"^ = 4 X 10^"^, which corresponds fairly well to 
the dynamic range of about 2.5 decades of the measure- 
ment. The same argument can be applied to the accuracy 
of the baseline in the power spectrum determined with 
the WTSA and presented in Fig. 4. In this case, how- 
ever, the baseline is generated by the shot noise, which 
is correlated on a time of the order of the inverse of the 
sampling frequency uj^. The number N of independent 
samples corresponds to the total number of samples ac- 
cumulated at ujQ. For the results presented in Fig. 4 
N = 10^ and the accuracy of the baseline is of the order 
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of l/iV"'^/^ = 3 X 10 ^. By eliminating the contribution 



of afterpulsing to the spectrum as described in Sect. IV 



the dynamic range could potentially be extended about 
five decades above the intrinsic contrast of the measure- 
ment, which in this case amounts to about 10'^. It has to 
be remarked that in this case there's no upper limit for 
the contrast, therefore a dynamic range as high as 10^° 
could potentially be obtained. 



B. Model system 

The ultimate task of the WTSA is the characteriza- 
tion of signals characterized by a wide distribution of 
timescales. 

To assess its performances we have chosen a sim- 
ple model system made up by strongly diluted diffus- 
ing brownian particles illuminated by a wide coUimated 
laser beam. This system gives rise to a nice power-law 
power spectrum with the exponent —3/2. A detailed 
derivation of the predicted power spectrum is presented 
in Appendix. The sample is a colloidal suspension of 
polystyrene latex spheres of 10/im diameter, suspended 
in water at a concentration of roughly 10'*particles/cm'^. 
The solvent is a mixture of equal volumes of water and 
deuterated water, so that its density is matched to that 
of the colloid within 0.1%. The sample is illuminated by 
a gaussian laser beam with a diameter at 1/e^ of about 
1mm. The scattering volume is delimited by two slits 
perpendicular to the beam. Therefore the scattering vol- 
ume is determined in one direction by the hard edge of 
the slits and in the other direction by the gaussian profile 
of the beam. The scattered light is collected at 90° with 
respect to the main beam. 

Due to the strongly diluted sample and to the wide 
beam illumination, fluctuations in the intensity of the 
scattered light are mostly determined by changes in the 
number of particles in the scattering volume. By using 
a wide laser beam, the coherence areas onto the photo- 
catode of the PMT are very small. Therefore, the pho- 
tocatode collects the averaged contribution coming from 
many coherence areas. In this way the component of 
the correlation function due to the diffusion of particles 
has a vanishingly small contrast, when compared to the 
averaged intensity measured by the photodetector. The 
signal onto the PMT thereby represents the superposi- 
tion of the intensities scattered by single particles. As it 
is shown in the Appendix, intensity fluctuations spanning 
a wide frequency range are excited when particles cross 
the hard edge at the boundary of the scattering volume. 

Figure 5 shows the correlation function of the fluctu- 
ations in the scattered intensity. The correlation func- 
tion shows an almost featureless structure and a narrow 
dynamic range. The correlation function is almost flat 
within about six decades in lag time, due to the long 
tails associated with the slowest modes. The peak at 
small lags still represents the contribution of afterpulses. 
The correlation function eventually decays at times of 
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FIG. 5: Correlation function of the intensity fluctuations in 
the model system, measured with ALV equipement. 



the order of 60s. However, this decay is an artifact due 
to the finite measurement time. In fact, correlation anal- 
ysis assumes that the correlation function decays to zero 
at large delays, once its baseline has been subtracted. 
The baseline is usually evaluated from the value of the 
correlation function at very large delays. However, for 
time-scale invariant processes, the baseline still contains 
the contribution of the slowest modes. Therefore, an in- 
crease of the acquisition time involves a decay at a larger 
time. As a rule of thumb, experimental data for the cor- 
relation function are reliable up to delays roughly 100 
times smaller than the acquisition time. This is a feature 
of correlation techniques which strongly limits their abil- 
ity in characterizing processes where an upper time scale 
is absent. 

Figure 6 shows the power spectrum of the fluctuations 
measured by means of the WTSA. The crosses represent 
the raw data, the horizontal line the background level 
due to shot noise and the circles the subtracted data rep- 
resenting the power spectrum. The most prominent fea- 
ture of the results in Fig. 6 is the impressive dynamic 
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FIG. 6: Power spectrum of the intensity fluctuations in the 
model system, measured with WTSA. Crosses represent the 
raw data. The dashed line is the baseline due to the shot noise. 
Circles represent the power spectrum obtained by subtracting 
the baseline from the raw data. The continuous line represents 
a power-law with exponent —3/2 for comparison. 

range spanning about 8 decades. The measured power 
spectrum decays approximately as a power law with an 
exponent of the order of —1.2. The straight line repre- 
sents a power law with exponent —3/2 for comparison. 
At the larger frequency, the contribution of afterpulses is 
apparent. 

Deviations from the power law behavior at intermedi- 
ate frequencies can be attributed to a residual contribu- 
tion of the diffusion of particles to the correlation func- 
tion. Such a contribution cannot be completely elimi- 
nated, even by collecting a large number of coherence ar- 
eas onto the multiplier. As outlined above, the diffusion 
of particles gives rise to fluctuations in the scattered in- 
tensity at 90°, characterized by lorentzian spectrum with 
a linewidth of the order of T = 20Hz. Deviations form 
the power-law behavior occur at frequencies roughly cen- 
tered around F. However, a more refined analysis of the 
data, including the evaluation of contributions due to dif- 



fusion and due to number fluctuations to the spectrum, 
is beyond the aim of this paper. 

The wavelet transform analysis is particularly effective 
for processes where the power law exponent a is larger 
than 1, like the model system just described. In fact, 
correlation techniques completely fail when applied to 
such systems. For exponents smaller than one the spec- 
tral behavior can be characterized either by measuring 
the power spectrum S (uj) cx 1/uj" or by measuring the 
correlation function, which also exhibits a power-law be- 
havior C (t) oc t^^". In the limit a 1 the correlation 
function becomes flat. For a > 1 the correlation func- 
tion is still flat and time independent. Naively speaking, 
this is due to the strong divergence of the power spec- 
trum at small frequency which behaves as a delta shaped 
spectrum. The diverging energy of the small frequency 
modes gets uniformly distributed along the lag time axis 
of the correlation function, which becomes flat. More 
rigorously, the correlation function can be obtained by 
Fourier transforming the power-law spectrum, regular- 
ized to avoid its divergence at w = 0: 

C (t) (X J lini — j- — — exp {~iujt)duj. (6) 

For a > 1, the spectrum has an integrable power law tail 
for w — > -l-oo. It approximates a delta function S (w) as 
e ^ 0. Therefore, for a > 1, C {t) is a constant: a cor- 
relation function like that shown in Fig. 3 is completely 
featureless and useless. 

An important application of the WTSA is the time de- 
pendent analisys of transient behavior during processes 
such as gelation ^ and colloidal aggregation |^ . In 
general, during such transients the highest frequencies 
are excited at first, and lower frequencies get gradually 
excited as time goes by. This introduces a characteristic 
lag time Tc marking the boundary between excited and 
non-excited modes. The proper characterization of this 
processes involves the measurement of the power spec- 
trum (or correlation function) from data sampled along 
a duration time T larger then Tc, so to be able to resolve 
Tc. In principle a longer T allows to measure the power 
spectrum or correlation function with a higher accuracy, 
due to the increase of the number of independent sample 
processed. However the lag time Tc changes during the 
transient. Therefore, a compromise has to be reached 
between the need to keep T not much larger than Tc and 
the desire to measure the spectrum with a high accuracy. 

As far as the WTSA is concerned, the same argument 
used to estimate the accuracy of the baseline can be used 
to get an order of magnitude estimate of the number of 
samplings needed to get a chosen accuracy. Suppose that 
we want to obtain the power spectrum at a frequency ujk 
with an accuracy better than a. As it was shown in sec- 
tion H, the power spectrum is calculated by averaging 

the squared wavelet coefficients: S {ojk) X^n I'-nl ■ 
recall here that the coefficients represent the wavelet 
coefficients evaluated at subsequent times n. Although 
this argument is not rigorous, for a generic random sig- 
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nal, it is reasonable to assume that the coefficients eval- 
uated at different times are independent. Therefore the 

1 11 

percentual error on S [iOk) scales as l/A^^ , where Nk 
is the number of samplings of the signal accumulated at 
the frequency Wfc. In this way, a 1% in 5* (w^) is achieved 
by sampling the signal up to frequencies of the order of 
lOOOOwfc. However, an order of magnitude estimate of the 
power spectrum can be obtained from a few samplings at 

Therefore, a very important feature of the wavelet 
analysis, when compared to correlation techniques, is its 
ability to characterize the power spectrum at small fre- 
quencies, without introducing artifacts related to the fi- 
nite acquisition time. This can be appreciated from Fig. 
6, where the smallest experimental angular frequency is 
of the order of 10~^Hz, roughly corresponding to the ac- 
quisition time of 2700s. The wavelet analysis only yields 
an order of magnitude estimate of the amplitude of this 
mode, due to the absence of a significative statistical 
sample. However, the wide dynamic range involved in 
the power spectrum makes it almost insensitive to back- 
ground subtraction at small frequencies, so that no char- 
acteristic times associate with the duration of the experi- 
ment are introduced in the subtraction. Therefore, when 
compared with correlation techniques, wavelet analysis 
allows a dramatic increase of the investigated frequency 
range, this increase amounting to almost two decades. 
This feature makes the WTSA the ideal tool to investi- 
gate slow processes, non-stationary processes, time-scale 
invariant processes and processes characterized by a wide 
distribution of time-scales in general. 
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scattered by each particle is proportional to the intensity 
/ (x) of the main beam at the position of the particle: 



It(X ^/[f„(i)]. 



(Al) 



n=0 



Roughly speaking. It is proportional to the number of 
particles inside the main beam, weighted by the intensity 
in that point. The temporal correlation function is thus: 



JV-l 

Cl {t) ^ {Ut+r) ^ J2 {I[^n{t)]I[£m{t + T)]) , 
n,m— 

(A2) 

where the brackets represent the average over t. By as- 
suming that the particles are non-interacting, their posi- 
tions are not correlated: 

Ci{t)(xN{N-1)P + N{I [x (t)] I[x{t + r)]) , (A3) 

where the bar represents the average over the volume 
V of the whole cell. From this equation, we derive the 
relative mean square value of the fluctuations of It : 



(If) - {Itf Ci (0) - Ci (+oo) 1 P - P 



(It 



Ci i+oo) 



N P 



(A4) 



Roughly speaking, (^P — P) /P is the ratio V/Vb, where 
Vb is the volume inside the main beam. This statement 
can be easily proved for an hypotetic beam with constant 
intensity inside Vb, vanishing outside. By calling Nb the 
number of particles inside the beam: 
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(A5) 
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APPENDIX A: POWER SPECTRUM OF THE 
NUMBER FLUCTUATIONS OF BROWNIAN 
PARTICLES IN A GIVEN VOLUME. 

In this appendix we will derive the power spectrum 
determined by the number fluctuations of brownian par- 
ticles diffusing within a scattering volume delimited by 
sharp edges on two sides and by the gaussian profile of 
a beam on the two perpendicular sides. We will assume 
that the sample is illuminated by a wide coUimated beam, 
so that spatial coherence effects can be neglected. Some 
of the arguments we use in the following are also used in 
the derivation of fluorescence correlation spectroscopy. 

Given N particles, at position Xn, the measured inten- 
sity /( at time t is the sum of the intensities scattered by 
each particle, since light is collected over several speckles, 
and interference plays no role. Moreover, the intensity 



This equation is analogous to the relation between num- 
ber fluctuation and mean number for Poisson distribu- 
tion. In order for the fluctuations to be strong, only a 
few particles must be inside the beam simultaneously. 

The diffusion of each particle can be described as a 
random walk; P^. {y — x), the probability that a particle 
in X goes in y in a time r, is thus: 



Pr (Ax) 



Ax" 
' 20-2 (r) 



(A6) 



where a (t) = y^D \t\ and D is the diffusion coefficient. 
By using thi s pr obability distribution, we can explicitly 
rewrite Eq. (A3): 



C/ (r) cx c -t- J I (x) I (y) {y - x) dxdy, ( A7) 

where c is a constant. The convolution operator is diag- 
onalized by Fourier transform: 



\Iiq)\'Priq)dq, 



(A8) 
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where Pr (q) and / (q) are the Fourier transforms in x of 
Pr {x) and / {x) : 



Pr{q) = e 2" ^ \ 



By Fourier transforming Eq. (A8) in r: 



(A9) 



Now, we impose that cu ^ DQj^ and w ^ DQ\. This 
means that w corresponds to times much shorter than 
the time needed by a particle to travel across the main 
beam. Under this condition, we can neglect qx and qy at 
the denominator, since the gaussian part in the integrand 
imposes |(?a;| ^ Qr and |(7y| ^ Qr', then, we integrate over 
qx and qy-. 



\I{q)\' P^{q)Aq, 



(AlO) 



where Puj {q) is the Fourier transform in r of P^ (q) as 
defined by Eq. @: 



(All) 



Now, we explicitly give the intensity distribution / (x) 
within the scattering volume: 



167r/g f sin^ jqJQL) D {ql + 2Q|) ^ 

(^) = ^2 ^2 / ^ ^rT2 ,.,2_Ln2^4 



(A15) 

The sin^ {qz/Qh) has fast oscillations, and averages to 
1/2 for q:, > Ql- 



-^^M±^d... (A16) 



D^qi 



I {x, y, z) = hX[-L/2X/2\ {z) e 



2i?2 , 



(A12) 



where L is the length of the portion of the beam we 
observe, R is the radius of the main beam at e"^/^ and 
Jo is the intensity at the center of the main beam. The 
Fourier transform is: 



J, ^ 47r/o sm[qzlQL) 

I{qx.qy.qz) = 77752- „ m ^ 

QlWr Qz/WL 



ql + ql 

2Ql 



(A13) 



where Ql = '^l L & Qr = 1/R- B y ins erting th e exp licit 
expressions of Eq. ( |A9| ) and Eq. ( A13D in Eq. ( [A1C| ): 



16^/g /■ sin^ jqJQL) 
''^^^^QWrJ (qz/Q.f ' 

u:^ + D^ (q^x + q'y + q'z) 



ql + ql 
Qr X (A14) 

jdqxdqydqz- 



The leading term in the integral is the coming from 
the sine function. This means that the behaviour of the 
measured intensity is mainly due to particles crossing the 
sharp edges of the slits, moving along the main beam, 
while the movements inside and outside the beam can be 
neglected. 

Since Qr and Ql are of the same order of magnitude: 



(A17) 



The integral converges; by substituting qz = xy^: 



(A18) 
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